package ua.lviv.franko.solvers.stressanalysis.twoDimention;

import java.util.ArrayList;
import java.util.Iterator;

import org.jzy3d.maths.Range;

import ua.lviv.franko.integration.ConstantFunction2D;
import ua.lviv.franko.integration.Gauss2D;
import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.matrix.RowType;
import ua.lviv.franko.matrix.hash.lu.HashMatrixLUStyle;
import ua.lviv.franko.matrix.hash.lu.LUFactorization;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.CubicSpline;
import ua.lviv.franko.mlsFunctions.twoDimension.BaseFunction;
import ua.lviv.franko.mlsFunctions.twoDimension.Particle;
import ua.lviv.franko.mlsFunctions.twoDimension.Values;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.IWeightFunction2D;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.SquareSupport;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.EquationParameters;

public class EFGOptimizedSolver {

	private static final double			EPS					= 0.000000000001;
	// general data
	protected int						quadraturePointsNumber;
	protected Gauss2D					gau;
	protected EquationParameters		parameters;
	protected HashMatrixLUStyle				mainMatrix;
	protected double[]					b;
	private static final int			NUMBER_OF_UNKNOWNS	= 2;

	// // data for background integration;
	// protected ArrayList<Rectangle> recList;
	// // protected ArrayList<Node> nodeList;
	// // protected ArrayList<NodePair> edges;
	// data for base functions
	protected ArrayList<Particle>		particles;
	protected ArrayList<BaseFunction>	baseFunctions;

	public EFGOptimizedSolver(int quadraturePointsNumber, ArrayList<Particle> particles, EquationParameters parameters, int PolinomialDegree, IWeightFunction2D weight) {
		// initialization
		this.quadraturePointsNumber = quadraturePointsNumber;
		this.gau = new Gauss2D(quadraturePointsNumber);
		this.particles = particles;
		this.parameters = parameters;

		// creation of base functions;
		this.createBaseFunctions(PolinomialDegree, weight);
		// continuation of initialization

	}

	public EFGResult2D calculate() {
		// for (int i = 0; i < b.length; i++) {
		// for (int j = 0; j < b.length; j++) {
		// System.out.print(mainMatrix.getValue(i, j)+"\t");
		// }
		// System.out.println();
		// }
		// for (int j = 0; j < b.length; j++) {
		// System.out.println(b[j]);
		// }
		for (int i = 0; i < baseFunctions.size(); i++) {
			BaseFunction bFunc = baseFunctions.get(i);
			// if ((Math.abs(Math.abs(bFunc.getXi().getXi()) - 0.5) < 0.000001)
			// || (Math.abs(Math.abs(bFunc.getXi().getYi()) - 0.5) < 0.000001))
			// {
			// // if ((Math.abs(Math.abs(bFunc.getXi().getYi()) - 0.5) <
			// // 0.0000001)) {
			// // for (int j = 0; j < mainMatrix.getN(); j++) {
			// // mainMatrix.setValue(i, j, 0.0);
			// // }
			// mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS, i
			// * NUMBER_OF_UNKNOWNS, 1000.0);
			// mainMatrix.setValue(i * NUMBER_OF_UNKNOWNS + 1, i
			// * NUMBER_OF_UNKNOWNS + 1, 1000.0);
			// b[i * NUMBER_OF_UNKNOWNS] = 0.0;
			// b[i * NUMBER_OF_UNKNOWNS] = 0.0;
			// }
			// upper
			double val = 0.0;
			if ((Math.abs(bFunc.getXi().getXi() - 0.5) < EPS)) {
				setBC(i, val);
				// continue;
			}
			// left
			val = 0.0;
			if ((Math.abs(bFunc.getXi().getXi() + 0.5) < EPS)) {
				setBC(i, val);
				// continue;
			}

		}
		LUFactorization lu = new LUFactorization(mainMatrix);
		System.out.println();
		double[] res = lu.calculate(b);
		// for (int i = 0; i < res.length; i++) {
		// System.out.println(res[i]);
		// }
		return new EFGResult2D(res, baseFunctions, particles);
		// return new EFGResult2D(res, baseFunctions);
	}

	public void generateMatrix() {
		mainMatrix = new HashMatrixLUStyle(RowType.ARRAY_TYPE, particles.size() * NUMBER_OF_UNKNOWNS);
		b = new double[particles.size() * NUMBER_OF_UNKNOWNS];

		for (int i = 0; i < b.length / NUMBER_OF_UNKNOWNS; i++) {
			// System.out.println(i);
			BaseFunction baseFunction = baseFunctions.get(i);
			double[][] vals = this.gaussQuadratureIntegration(2, Math.max(baseFunction.getXi().getXi() - baseFunction.getXi().getRho(), -0.5),
					Math.min(baseFunction.getXi().getXi() + baseFunction.getXi().getRho(), 0.5), Math.max(baseFunction.getXi().getYi() - baseFunction.getXi().getRho(), -0.5),
					Math.min(baseFunction.getXi().getYi() + baseFunction.getXi().getRho(), 0.5), baseFunction);
			// double[][] vals = this.gaussQuadratureIntegration(-0.5, 0.5,
			// -0.5, 0.5, baseFunction);
			int pos1 = particles.indexOf(baseFunction.getXi());
			for (int j = 0; j < vals.length; j++) {
				int pos2 = particles.indexOf(baseFunction.getNodeList().get(j));
				mainMatrix.setValue(pos1 * NUMBER_OF_UNKNOWNS, pos2 * NUMBER_OF_UNKNOWNS, vals[j][0]);
				mainMatrix.setValue(pos1 * NUMBER_OF_UNKNOWNS, pos2 * NUMBER_OF_UNKNOWNS + 1, vals[j][1]);
				mainMatrix.setValue(pos1 * NUMBER_OF_UNKNOWNS + 1, pos2 * NUMBER_OF_UNKNOWNS, vals[j][2]);
				mainMatrix.setValue(pos1 * NUMBER_OF_UNKNOWNS + 1, pos2 * NUMBER_OF_UNKNOWNS + 1, vals[j][3]);
			}
		}

		for (int i = 0; i < b.length; i++) {
			for (int j = 0; j < b.length; j++) {
				System.out.print(mainMatrix.getValue(i, j) + "\t");
			}
			System.out.println();
		}
		for (int j = 0; j < b.length; j++) {
			System.out.println(b[j]);
		}

	}

	private void createBaseFunctions(int PolinomialDegree, IWeightFunction2D weight) {
		this.baseFunctions = new ArrayList<BaseFunction>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle particle = iterator.next();
			BaseFunction bf = new BaseFunction(PolinomialDegree, particle, weight);
			bf.setNodeList(this.getParticleSupport(particle, weight));
			this.baseFunctions.add(bf);
			// bf.display();
		}
	}

	private double[][] gaussQuadratureIntegration(int n, double ax, double bx, double ay, double by, BaseFunction baseFunction) {
		double hx = (bx - ax) / n;
		double hy = (by - ay) / n;

		double[][] res = new double[baseFunction.getNodeList().size()][NUMBER_OF_UNKNOWNS * NUMBER_OF_UNKNOWNS];

		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				double[][] tmp = gaussQuadratureIntegration(ax + i * hx, ax + (i + 1) * hx, ay + j * hy, ay + (j + 1) * hy, baseFunction);
				for (int k = 0; k < tmp.length; k++) {
					for (int q = 0; q < tmp[0].length; q++) {
						res[k][q] += tmp[k][q];
					}
				}
			}
		}

		return res;
	}

	private double[][] gaussQuadratureIntegration(double ax, double bx, double ay, double by, BaseFunction baseFunction) {
		// first column - baseFunction*v; second column - v*baseFunction
		double[][] result = new double[baseFunction.getNodeList().size()][NUMBER_OF_UNKNOWNS * NUMBER_OF_UNKNOWNS];
		// its a position of main node - node of base function
		int pos = baseFunction.getNodeList().indexOf(baseFunction.getXi());

		double f1 = 0.0;
		double f2 = 0.0;

		for (int j = 0; j < quadraturePointsNumber; j++) {
			double temp1 = 0.0;
			double temp2 = 0.0;
			for (int i = 0; i < quadraturePointsNumber; i++) {
				double x = (ax + bx) / 2 + (bx - ax) / 2 * gau.getX(i);
				double y = (ay + by) / 2 + (by - ay) / 2 * gau.getX(j);
				Values[] values = baseFunction.calculateAllFunctions(x, y);
				double k1 = parameters.getK1();
				double k2 = parameters.getK2();
				double k3 = parameters.getK3();

				for (int k = 0; k < values.length; k++) {
					result[k][0] += gau.getW(j) * gau.getW(i) * (values[pos].valX * values[k].valX * k1 + values[pos].valY * values[k].valY * k3);
					result[k][1] += gau.getW(j) * gau.getW(i) * (values[pos].valY * values[k].valX * k2 + values[pos].valX * values[k].valY * k3);
					result[k][2] += gau.getW(j) * gau.getW(i) * (values[pos].valY * values[k].valX * k2 + values[pos].valX * values[k].valY * k3);
					result[k][3] += gau.getW(j) * gau.getW(i) * (values[pos].valX * values[k].valX * k3 + values[pos].valY * values[k].valY * k1);

				}

				temp1 += gau.getW(i) * values[pos].val * parameters.getF1().calculate(x, y);
				temp2 += gau.getW(i) * values[pos].val * parameters.getF2().calculate(x, y);
			}
			f1 += temp1 * gau.getW(j);
			f2 += temp2 * gau.getW(j);
		}
		f1 *= (bx - ax) * (by - ay) / 4.0;
		f2 *= (bx - ax) * (by - ay) / 4.0;
		for (int i = 0; i < result.length; i++) {
			result[i][0] *= (bx - ax) * (by - ay) / 4.0;
			result[i][1] *= (bx - ax) * (by - ay) / 4.0;
			result[i][2] *= (bx - ax) * (by - ay) / 4.0;
			result[i][3] *= (bx - ax) * (by - ay) / 4.0;
		}
		b[particles.indexOf(baseFunction.getXi()) * NUMBER_OF_UNKNOWNS] += f1;
		b[particles.indexOf(baseFunction.getXi()) * NUMBER_OF_UNKNOWNS + 1] += f2;
		return result;
	}

	private ArrayList<Particle> getParticleSupport(Particle p, IWeightFunction2D weight) {
		ArrayList<Particle> list = new ArrayList<Particle>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle particle = (Particle) iterator.next();
			if (weight.intersectedSupport(p, particle)) {
				list.add(particle);
			}
		}
		return list;
	}

	private void setBC(int pos, double val) {
		for (int i = 0; i < mainMatrix.getN(); i++) {
			mainMatrix.setValue(2 * pos, i, 0.0);
			mainMatrix.setValue(2 * pos + 1, i, 0.0);
		}
		mainMatrix.setValue(pos * 2, pos * 2, 10000.0);
		mainMatrix.setValue(pos * 2 + 1, pos * 2 + 1, 10000.0);
		b[2 * pos] = val * 10000;
		b[2 * pos + 1] = val * 10000;
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		ArrayList<Particle> nodes = new ArrayList<Particle>();
		// параметри області:
		double ax = -0.5, bx = 0.5, ay = -0.5, by = 0.5;

		int n = 10, m = 10;

		double hx = (bx - ax) / (double) n, hy = (by - ay) / (double) m;

		// вузли:
		double rho = 1.5 * hx;
		for (int i = 0; i <= n; i++) {
			for (int j = 0; j <= m; j++) {
				Particle p = new Particle(ax + i * hx, ay + j * hy, rho);
				nodes.add(p);
			}
		}

		// equation parameters:
		EquationParameters params = new EquationParameters();
		IFunction2D f1 = new ConstantFunction2D(1.0);
		IFunction2D f2 = new ConstantFunction2D(0.0);
		params.setF1(f1);
		params.setF2(f2);
		double e = 0.6;
		double v = 0.6;
		params.setK1(e / (1 - v * v));
		params.setK2(e * v / (1 - v * v));
		params.setK3(e / (1 - v * v) * ((1 - v) / 2.0));

		EFGOptimizedSolver solver = new EFGOptimizedSolver(5, nodes, params, 1, new SquareSupport(new CubicSpline()));

		solver.generateMatrix();
		solver.calculate().displayU(new Range(-0.5, 0.5), new Range(-0.5, 0.5));

	}
}
